home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / machserver / 1.098 / mach / sun4c.md / div.c < prev    next >
C/C++ Source or Header  |  1989-08-17  |  11KB  |  414 lines
  1. #ifndef lint
  2. static char     sccsid[] = "@(#)div.c 1.6 88/02/08 SMI";
  3. #endif
  4.  
  5. /*
  6.  * Copyright (c) 1988 by Sun Microsystems, Inc.
  7.  */
  8.  
  9. /*
  10.  * This whole file refuses to lint.
  11.  */
  12. #ifndef lint
  13. #include <sun4/fpu/fpu_simulator.h>
  14. #include <sun4/fpu/globals.h>
  15.  
  16. PRIVATE unsigned
  17. trial_quotient(remainder, divisor)
  18.     unsigned        remainder[], divisor[];
  19.  
  20. /*
  21.  * Given remainder[0,1,2,3,4] < 2**16, divisor[0,1,2,3] < 2**16, 2**15 <=
  22.  * divisor[0] < 2**16, computes a trial quotient q < 2**17, such that q =
  23.  * remainder/divisor, and revises remainder.
  24.  */
  25.  
  26. {
  27.     unsigned        q;
  28.     int             i, k;
  29.     unsigned        t, t1, t2;
  30.     q = ((remainder[0] << 16) + remainder[1]) / divisor[0];
  31.  
  32.     for (i = 3; i >= 0; i--) {
  33.         t = (q & 0xffff) * divisor[i];
  34.         t1 = t >> 16;
  35.         t2 = t & 0xffff;
  36.  
  37.         if (remainder[i + 1] < t2) {    /* ripple carry */
  38.             k = i;
  39.             while ((k >= 0) && (remainder[k] == 0))
  40.                 remainder[k--] = 0xffff;
  41.             remainder[k]--;
  42.         }
  43.         remainder[i + 1] -= t2;
  44.         remainder[i + 1] &= 0xffff;
  45.  
  46.         if ((i >= 1) && (remainder[i] < t1)) {    /* ripple carry */
  47.             k = i - 1;
  48.             while ((k >= 0) && (remainder[k] == 0))
  49.                 remainder[k--] = 0xffff;
  50.             remainder[k]--;
  51.         }
  52.         remainder[i] -= t1;
  53.         remainder[i] &= 0xffff;
  54.  
  55.         t = (q >> 16) * divisor[i];
  56.         if (t > 0) {
  57.             if ((i >= 1) && (remainder[i] < t)) {    /* ripple carry */
  58.                 k = i - 1;
  59.                 while ((k >= 0) && (remainder[k] == 0))
  60.                     remainder[k--] = 0xffff;
  61.                 remainder[k]--;
  62.             }
  63.             remainder[i] -= t;
  64.             remainder[i] &= 0xffff;
  65.         }
  66.     }
  67.  
  68.     while (remainder[0] > 0x8000) {    /* Sign reversal occurred - reduce
  69.                      * quotient by 1 - add divisor back
  70.                      * to remainder. */
  71.         q--;
  72.         for (i = 3; i >= 0; i--) {
  73.             if (remainder[i + 1] > (0xffff - divisor[i])) {    /* ripple carry */
  74.                 k = i;
  75.                 while ((k >= 0) && (remainder[k] == 0xffff))
  76.                     remainder[k--] = 0;
  77.                 if (k >= 0)
  78.                     remainder[k]++;
  79.             }
  80.             remainder[i + 1] += divisor[i];
  81.             remainder[i + 1] &= 0xffff;
  82.         }
  83.     }
  84.  
  85.     /* Make sure that remainder <= divisor. */
  86.     i = 0;
  87.     while (remainder[i + 1] == divisor[i])
  88.         i++;
  89.     if ((remainder[0] != 0) || (i > 3) || (remainder[i + 1] > divisor[i])) {    /* Remainder too big -
  90.                                              * increase quotient by
  91.                                              * 1 - subtract divisor
  92.                                              * back from remainder. */
  93.         q++;
  94.         for (i = 3; i >= 0; i--) {
  95.             if (remainder[i + 1] < divisor[i]) {    /* ripple carry */
  96.                 k = i;
  97.                 while ((k >= 0) && (remainder[k] == 0))
  98.                     remainder[k--] = 0xffff;
  99.                 if (k >= 0)
  100.                     remainder[k]--;
  101.             }
  102.             remainder[i + 1] -= divisor[i];
  103.             remainder[i + 1] &= 0xffff;
  104.         }
  105.         i = 0;
  106.         while (remainder[i] == divisor[i])
  107.             i++;
  108.     }
  109.     for (i = 0; i <= 3; i++)
  110.         remainder[i] = remainder[i + 1];
  111.     remainder[4] = 0;
  112.  
  113.     return q;
  114. }
  115.  
  116. void
  117. _fp_div(px, py, pz)
  118.     unpacked       *px, *py, *pz;
  119.  
  120. {
  121.     unsigned        quotient[5];    /* Quotient accumulator. */
  122.     unsigned        remainder[5];    /* Remainder accumulator. */
  123.     unsigned        divisor[5];
  124.  
  125.     *pz = *px;
  126.     pz->sign = px->sign ^ py->sign;
  127.  
  128.     if ((py->fpclass == fp_quiet) || (py->fpclass == fp_signaling)) {
  129.         *pz = *py;
  130.         return;
  131.     }
  132.     switch (px->fpclass) {
  133.     case fp_quiet:
  134.     case fp_signaling:
  135.         return;
  136.     case fp_zero:
  137.     case fp_infinity:
  138.         if (px->fpclass == py->fpclass) {    /* 0/0 or inf/inf */
  139.             fpu_error_nan(pz);
  140.             pz->fpclass = fp_quiet;
  141.         }
  142.         return;
  143.     case fp_normal:
  144.         switch (py->fpclass) {
  145.         case fp_zero:    /* number/0 */
  146.             fpu_set_exception(fp_division);
  147.             pz->fpclass = fp_infinity;
  148.             return;
  149.         case fp_infinity:    /* number/inf */
  150.             pz->fpclass = fp_zero;
  151.             return;
  152.         }
  153.     }
  154.  
  155.     /* Now x and y are both normal or subnormal. */
  156.  
  157.     pz->exponent = px->exponent - py->exponent;
  158.  
  159.     remainder[0] = px->significand[0] >> 16;
  160.     remainder[1] = px->significand[0] & 0xffff;
  161.     remainder[2] = px->significand[1] >> 16;
  162.     remainder[3] = px->significand[1] & 0xffff;
  163.     remainder[4] = 0;
  164.     divisor[0] = py->significand[0] >> 16;
  165.     divisor[1] = py->significand[0] & 0xffff;
  166.     divisor[2] = py->significand[1] >> 16;
  167.     divisor[3] = py->significand[1] & 0xffff;
  168.     divisor[4] = 0;
  169.  
  170.     quotient[0] = trial_quotient(remainder, divisor);
  171.     quotient[1] = trial_quotient(remainder, divisor);
  172.     pz->significand[0] = (quotient[0] << 16) + quotient[1];
  173.     quotient[2] = trial_quotient(remainder, divisor);
  174.     quotient[3] = trial_quotient(remainder, divisor);
  175.     pz->significand[1] = (quotient[2] << 16) + quotient[3];
  176.     if (quotient[2] > 0xffff)
  177.         pz->significand[0]++;
  178.  
  179.     if (remainder[0] >= 0x8000) {    /* remainder is large enough to force
  180.                      * round upward. */
  181.         pz->significand[2] = 0xc0000000;
  182.     } else {        /* Remainder < 1.0 <= divisor so compare
  183.                  * 2*remainder to divisor. */
  184.         int             i;
  185.         for (i = 0; i <= 3; i++) {
  186.             remainder[i] <<= 1;
  187.             if (remainder[i] >= 0x10000) {    /* Shift out. */
  188.                 remainder[i] &= 0xffff;
  189.                 remainder[i - 1] |= 1;
  190.             }
  191.         }
  192.         for (i = 0; i <= 3; i++) {    /* compare divisor to
  193.                          * 2*remainder */
  194.             if (divisor[i] > remainder[i])
  195.                 goto gt;
  196.             if (divisor[i] < remainder[i])
  197.                 goto lt;
  198.         }
  199.         /*
  200.          * Can't fall through to here - remainder can't be exactly
  201.          * half divisor.
  202.          */
  203.  
  204. gt:                /* divisor is greater than 2 * remainder */
  205.         pz->significand[2] = remainder[0] | remainder[1] | remainder[2] | remainder[3];
  206.         if (pz->significand[2] & 1)
  207.             pz->significand[2] |= 1;    /* Save lsb. */
  208.         pz->significand[2] >>= 1;    /* Right shift for leading 0
  209.                          * round bit; rest of word
  210.                          * will be zero if remainder
  211.                          * is exactly zero. */
  212.         goto round;
  213.  
  214. lt:                /* divisor is less than 2 * remainder. */
  215.         pz->significand[2] = 0xc0000000;    /* Therefore round
  216.                              * upward. */
  217.         goto round;
  218.     }
  219.  
  220. round:
  221.  
  222.     if (quotient[0] > 0xffff) {    /* 1 <= quotient < 2 */
  223.         fpu_rightshift(pz, 1);
  224.         pz->significand[0] |= 0x80000000;
  225.     } else
  226.         /* 0.5 <= quotient < 1 */
  227.         pz->exponent--;
  228.  
  229. }
  230.  
  231. #define ADDBIT(X,I) /* Adds  2**-i to x, an array of 16-bit chunks. */ \
  232.     (X)[((I)-1) / 16] += 0x8000 >> ((I-1) % 16)
  233.  
  234. #define SUBBIT(X,I) /* Adds -2**-i to x, an array of 16-bit chunks. */ \
  235.     (X)[((I)-1) / 16] -= 0x8000 >> ((I-1) % 16)
  236.  
  237.  
  238. PRIVATE
  239. dosqrt(px, pz)
  240.     unpacked       *px, *pz;
  241. {
  242.     int             i, j, negative, borrow;
  243.     unsigned        quotient[5];    /* Quotient accumulator. */
  244.     unsigned        remainder[5];    /* Remainder accumulator. */
  245.  
  246.     for (i = 0; i <= 4; i++) {
  247.         quotient[i] = 0;
  248.     }
  249.     remainder[0] = (px->significand[0] - 0x40000000) >> 16;
  250.     remainder[1] = px->significand[0] & 0xffff;
  251.     remainder[2] = px->significand[1] >> 16;
  252.     remainder[3] = px->significand[1] & 0xffff;
  253.     remainder[4] = px->significand[2] >> 16;
  254.  
  255.     negative = 0;
  256.     for (i = 1; i <= 65; i++) {    /* Generate one more bit of quotient. */
  257.         for (j = 0; j <= 4; j++) {    /* Double remainder. */
  258.             remainder[j] <<= 1;
  259.         }
  260.         if (negative) {    /* Remainder was negative - add quotient. */
  261.             for (j = 0; j <= 4; j++) {
  262.                 remainder[j] -= quotient[j];
  263.             }
  264.             SUBBIT(remainder, i + 1);
  265.             SUBBIT(remainder, i + 2);
  266.             for (j = 4; j >= 1; j--) {    /* Force
  267.                              * borrows/carries. */
  268.                 while (remainder[j] > 0x20000) {
  269.                     remainder[j] += 0x10000;
  270.                     remainder[j - 1] -= 1;
  271.                 }
  272.                 while (remainder[j] > 0xffff) {
  273.                     remainder[j] -= 0x10000;
  274.                     remainder[j - 1] += 1;
  275.                 }
  276.             }
  277.             borrow = 0;
  278.             while (remainder[0] > 0xffff) {
  279.                 remainder[0] += 0x10000;
  280.                 borrow++;
  281.             }
  282.             if ((remainder[0] | remainder[1] | remainder[2] | remainder[3] | remainder[4]) == 0) {    /* Remainder is exactly
  283.                                                          * zero. */
  284.                 negative = 0;
  285.             } else if (borrow) {    /* Remainder has changed
  286.                          * sign. */
  287.                 negative = 0;
  288.                 j = 4;
  289.                 while (remainder[j] == 0)
  290.                     j--;
  291.                 remainder[j] = 0x10000 - remainder[j];
  292.                 for (j--; j >= 0; j--) {    /* Complement remainder. */
  293.                     remainder[j] = 0xffff - remainder[j];
  294.                 }
  295.             }
  296.         } else {    /* Remainder is positive - subtract quotient. */
  297.             ADDBIT(quotient, i);    /* Adjust quotient. */
  298.             for (j = 0; j <= 4; j++) {
  299.                 remainder[j] -= quotient[j];
  300.             }
  301.             SUBBIT(remainder, i + 2);
  302.             for (j = 4; j >= 1; j--) {    /* Force
  303.                              * borrows/carries. */
  304.                 while (remainder[j] > 0x20000) {
  305.                     remainder[j] += 0x10000;
  306.                     remainder[j - 1] -= 1;
  307.                 }
  308.                 while (remainder[j] > 0xffff) {
  309.                     remainder[j] -= 0x10000;
  310.                     remainder[j - 1] += 1;
  311.                 }
  312.             }
  313.             borrow = 0;
  314.             while (remainder[0] > 0xffff) {
  315.                 remainder[0] += 0x10000;
  316.                 borrow++;
  317.             }
  318.             if (borrow) {    /* Remainder has changed sign. */
  319.                 negative = 1;
  320.                 j = 4;
  321.                 while (remainder[j] == 0)
  322.                     j--;
  323.                 remainder[j] = 0x10000 - remainder[j];
  324.                 for (j--; j >= 0; j--) {    /* Complement remainder. */
  325.                     remainder[j] = 0xffff - remainder[j];
  326.                 }
  327.             }
  328.         }
  329.     }
  330.     /* Quotient so far is exact iff remainder+quotient+2**-67 <= 0. */
  331.     if (negative) {        /* Maybe exact. */
  332.         for (j = 0; j <= 4; j++) {
  333.             remainder[j] -= quotient[j];
  334.         }
  335.         SUBBIT(remainder, 67);
  336.         for (j = 4; j >= 1; j--) {    /* Force borrows. */
  337.             while (remainder[j] > 0xffff) {
  338.                 remainder[j] += 0x10000;
  339.                 remainder[j - 1] -= 1;
  340.             }
  341.         }
  342.         if ((remainder[1] | remainder[2] | remainder[3] | remainder[4]) == 0) {    /* Remainder is exactly
  343.                                              * zero. */
  344.             negative = 0;
  345.         } else if (remainder[0] > 0xffff) {    /* Remainder has changed
  346.                              * sign. */
  347.             negative = 0;
  348.             j = 4;
  349.             while (remainder[j] == 0)
  350.                 j--;
  351.             remainder[j] = -remainder[j];
  352.             for (; j >= 0; j--) {    /* Complement remainder. */
  353.                 remainder[j] = 0xffffffff - remainder[j];
  354.             }
  355.         }
  356.     }
  357.     if (!negative &&
  358.         ((remainder[0] | remainder[1] | remainder[2] | remainder[3] | remainder[4]) != 0)) {    /* Remainder is neither
  359.                                                      * negative or zero and
  360.                                                      * therefore inexact. */
  361.         quotient[4] |= 1;
  362.     }
  363.     pz->significand[0] = ((quotient[0]) << 16) | quotient[1];
  364.     pz->significand[1] = ((quotient[2]) << 16) | quotient[3];
  365.     pz->significand[2] = ((quotient[4]) << 16);
  366. }
  367.  
  368. void
  369. _fp_sqrt(px, pz)
  370.     unpacked       *px, *pz;
  371.  
  372. {                /* *pz gets sqrt(*px) */
  373.  
  374.  
  375.     *pz = *px;
  376.     switch (px->fpclass) {
  377.     case fp_quiet:
  378.     case fp_signaling:
  379.     case fp_zero:
  380.         return;
  381.     case fp_infinity:
  382.         if (px->sign == 1) {    /* sqrt(-inf) */
  383.             fpu_error_nan(pz);
  384.             pz->fpclass = fp_quiet;
  385.         }
  386.         return;
  387.     case fp_normal:
  388.         if (px->sign == 1) {    /* sqrt(-norm) */
  389.             fpu_error_nan(pz);
  390.             pz->fpclass = fp_quiet;
  391.             return;
  392.         }
  393.     }
  394.  
  395.     /* Now x is normal. */
  396.  
  397.     if (px->exponent & 1) {    /* sqrt(1.f * 2**odd) = sqrt (0.5f) *
  398.                  * 2**(odd+1)/2 */
  399.         pz->exponent = (px->exponent - 1) / 2;
  400.     } else {        /* sqrt(1.f * 2**even) = sqrt (1.f) *
  401.                  * 2**(even)/2 */
  402.         fpu_rightshift(px, 1);    /* denormalize by 1 */
  403.         pz->exponent = px->exponent / 2;
  404.     }
  405.     /*
  406.      * Exponent is now correct unless sqrt rounds up to 2.0 - which can
  407.      * happen in round to positive infinity mode.
  408.      */
  409.  
  410.     dosqrt(px, pz);
  411.  
  412. }
  413. #endif /* lint */
  414.